#coding=utf-8
import argparse,os
import torch
from torch.autograd import Variable
import numpy as np
import time, math, glob
import scipy.io as sio
import matplotlib.pyplot as plt
import h5py

parser = argparse.ArgumentParser(description="PyTorch LapSRN Eval")
parser.add_argument("--cuda", default=True,action="store_true", help="use cuda?")
parser.add_argument("--model", default="/home/wks/users/lxj/MRI/pytorch-LapSRN-master/lapsrn_unet4.3_DC_0.0001_50/", type=str, help="model path")
parser.add_argument("--dataset", default="Set5", type=str, help="dataset name, Default: Set5")
parser.add_argument("--scale", default=4, type=int, help="scale factor, Default: 4")
parser.add_argument("--start_epoch",default=36,type=int)
parser.add_argument("--max_epoch",default=44,type=int)
parser.add_argument("--gpus", default="0", type=str, help="gpu ids (default: 0)")

def PSNR_3D(pred, gt):
    maxv = 1.0
    res = pred - gt
    # print(sum(res.flatten()))
    rmse = math.sqrt(np.mean(res ** 2))
    # print(f'rmse = {rmse}')
    if rmse == 0:
        return 100, 0
    psnr = 20 * math.log10(maxv / rmse)
    return psnr, rmse

def forward(t1,T2Dn,model):
    t1 = Variable(torch.from_numpy(t1).contiguous().float())
    t1= t1.view(1, -1, t1.shape[-2], t1.shape[-1])
    T2Dn = Variable(torch.from_numpy(T2Dn).contiguous().float())
    T2Dn = T2Dn.view(1, -1, T2Dn.shape[-2], T2Dn.shape[-1])
    
    mask = np.zeros((T2Dn.shape))
    d_i = range(148,190)
    mask[:,:,d_i] = 1
    mask_v = Variable(torch.from_numpy(mask).float())
    mask_v = mask_v.view(mask_v.shape[0],mask_v.shape[1],mask_v.shape[2],mask_v.shape[3],1)
    mask_v2 = torch.cat((mask_v,mask_v),-1)

    t1 = t1.cuda()
    T2Dn = T2Dn.cuda()
    mask_v2 = mask_v2.cuda()
#         x = x.to(self.device)

#     y = model(x) #网络应该返回两个结果
    y1,y2 = model(t1,T2Dn,mask_v2) #网络应该返回两个结果
    
    y1 = y1.cpu().data[0].numpy()  # .astype(np.float64) # maybe 1 channel
    y1 = y1.clip(0, 1.0)  # !!!!!!
    y2 = y2.cpu().data[0].numpy()  # .astype(np.float64) # maybe 1 channel
    y2 = y2.clip(0, 1.0)  # !!!!!!
    
#     print("y.cpu().data.numpy.shape:",y.cpu().data.numpy().shape)
#     y = y.cpu().data[0].numpy()  # .astype(np.float64) # maybe 1 channel
#     y = y.clip(0, 1.0)  # !!!!!!
    return y1,y2

def predict(model,imslices,T2Dn, slice_h, real_imsize):#slice_h=3
#     print("--------进入model_analy.tester.py---------------------")
    y_dim, z_dim = imslices.shape[2:]
        
    imsize = (imslices.shape[0] + slice_h - 1, y_dim, z_dim)

#     Y1, C = np.zeros(imsize, np.single),  np.zeros(imsize, np.single)
    Y1,Y2, C = np.zeros(imsize, np.single),np.zeros(imsize, np.single),  np.zeros(imsize, np.single)

    for i in range(imslices.shape[0]):
        t1 = imslices[i, :]
        t2Dn = T2Dn[i,:]

#         y= forward(x,model)

        y1,y2= forward(t1,t2Dn,model)
        
        C[i: (i + slice_h), :] += 1
        Y1[i: (i + slice_h), :] += y1
        Y2[i: (i + slice_h), :] += y2
    Y1 = Y1 / C
    Y1 = Y1[:, :real_imsize[1], :real_imsize[2]]
    Y2 = Y2 / C
    Y2 = Y2[:, :real_imsize[1], :real_imsize[2]]
    print("--------退出model_analy.tester.py---------------------")
    return Y1,Y2

opt = parser.parse_args()
cuda = opt.cuda

# if cuda and not torch.cuda.is_available():
#     raise Exception("No GPU found, please run without --cuda")
if cuda:                                                   
        print("=> use gpu id: '{}'".format(opt.gpus))            
        os.environ["CUDA_VISIBLE_DEVICES"] = opt.gpus
        if not torch.cuda.is_available():
                raise Exception("No GPU found or Wrong gpu id, please run without --cuda")


y_psnr1 = []
y_psnr2 = []
y_ssim = []
x = []
best_psnr1 = 0
best_psnr2 = 0
best_ssim = 0
psnr_i=ssim_i = 0
psnr_i2=ssim_i2 = 0

h5file = h5py.File("/home/wks/users/lxj/MRI/pytorch-LapSRN-master/data/S08.Sub(1).Slice3.T2-T1-T2Dn8.lapSRN2.Test..h5","r")
train_T1 = h5file["T1"]
T2Dn = h5file["T2Dn"]
label = h5file['label_T2']
label_T2Dn4 = h5file['label_Dn4']
T1 = h5file['t1_list']

for i in range(opt.start_epoch,opt.max_epoch+1):          #进入模型循环
    model_path = opt.model + "lapsrn_UNet_model_epoch_{}.pth".format(i)
    print(model_path)
#     model = torch.load(model_path, map_location=lambda storage, loc: storage)["model"]    
    model = torch.load(model_path)["model"]
    if cuda:
        model = model.cuda()        
    else:
        model = model.cpu()
    
    T2Dn4_HR,T2_HR = predict(model, train_T1,T2Dn, 3, label.shape)
    
    T2_HR[T1 == 0] = 0
    T2Dn4_HR[T1 == 0] = 0
    psnr, rmse = PSNR_3D(T2_HR, label)
    psnr2, rmse2 = PSNR_3D(T2Dn4_HR, label_T2Dn4)
    
    print("epoch={},psnr2={},psnr={}:".format(i,psnr2,psnr))
#     ssim = compare_ssim(T2_HR, label)
    y_psnr1.append(psnr)
    y_psnr2.append(psnr2)
#     y_ssim.append(ssim)

    x.append(i)
    if psnr > best_psnr1:
        best_psnr1 = psnr
        psnr_i = i
    if psnr2 > best_psnr2:
        best_psnr2 = psnr2
        psnr_i2 = i
y_psnr1 = np.array(y_psnr1)
y_psnr2 = np.array(y_psnr2)
data_file1 = opt.model + "PSNR1Test.npz"
data_file2 = opt.model + "PSNR2Test.npz"
np.savez(data_file1, psnr=y_psnr1)
np.savez(data_file2, psnr=y_psnr2)
print("Save Done")
# for image_name in image_list:
#     print("Processing ", image_name)
#     im_gt_y = sio.loadmat(image_name)['im_gt_y']
#     im_b_y = sio.loadmat(image_name)['im_b_y']
#     im_l_y = sio.loadmat(image_name)['im_l_y']
# 
#     im_gt_y = im_gt_y.astype(float)
#     im_b_y = im_b_y.astype(float)
#     im_l_y = im_l_y.astype(float)
# 
#     psnr_bicubic = PSNR(im_gt_y, im_b_y,shave_border=opt.scale)
#     avg_psnr_bicubic += psnr_bicubic
# 
#     im_input = im_l_y/255.
# 
#     im_input = Variable(torch.from_numpy(im_input).float()).view(1, -1, im_input.shape[0], im_input.shape[1])
# 
#     if cuda:
#         model = model.cuda()
#         im_input = im_input.cuda()
#     else:
#         model = model.cpu()
# 
#     start_time = time.time()
#     HR_2x, HR_4x = model(im_input)
#     elapsed_time = time.time() - start_time
#     avg_elapsed_time += elapsed_time
# 
#     HR_4x = HR_4x.cpu()
# 
#     im_h_y = HR_4x.data[0].numpy().astype(np.float32)
# 
#     im_h_y = im_h_y*255.
#     im_h_y[im_h_y<0] = 0
#     im_h_y[im_h_y>255.] = 255.
#     im_h_y = im_h_y[0,:,:]
# 
#     psnr_predicted = PSNR(im_gt_y, im_h_y,shave_border=opt.scale)
#     avg_psnr_predicted += psnr_predicted

# print("Scale=", opt.scale)
# print("Dataset=", opt.dataset)
# print("PSNR_predicted=", avg_psnr_predicted/len(image_list))
# print("PSNR_bicubic=", avg_psnr_bicubic/len(image_list))
# print("It takes average {}s for processing".format(avg_elapsed_time/len(image_list)))
print("max_epoch={},max_psnr2={}".format(psnr_i2,best_psnr2))
print("max_epoch={},max_psnr={}".format(psnr_i,best_psnr1))